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FINAL  REPORT 


MAGNETO-FLUID  DYNAMICS  CALCULATIONS  FOR  AERODYNAMICS 
AFOSR  GRANT  FA9550-04-1-0155 

Robert  W.  MacCormack 
Department  of  Aeronautics  and  Astronautics 
Stanford  University 

Abstract 

Magneto-Hydro-Dynamics  (MHD)  or  more  generally  Magneto-Fluid-Dynamics  (MFD)  offers  a 
potential  breakthrough  in  both  hypersonic  vehicle  design  and  propulsion.  Reductions  in  heat 
transfer  and  flow  control  using  magnetic  fields  can  be  important  for  enabling  a  hypersonic 
vehicle  to  pass  more  efficiently  and  safely  through  the  atmosphere.  Magnetic  and  electric  fields 
placed  within  the  propulsion  system  may  enable  the  extraction  of  electrical  energy  from  the 
ionized  flow  entering  the  engine,  while  simultaneously  slowing  the  fluid,  without  losses  in  total 
pressure  caused  by  shockwaves,  and  enhancing  complete  fuel  combustion.  The  extracted  energy 
can  be  returned  back  into  the  flow  after  combustion  for  further  flow  acceleration  and  engine 
thrust.  These  potential  benefits  may  or  may  not  be  realizable.  Realistic  aerodynamic  simulations, 
under  the  conditions  of  expected  low  electrical  conductivities  and  strong  magnetic  fields,  will  be 
required.  The  solution  of  the  complete  equations  governing  magneto-fluid  dynamics,  including 
magnetic  induction  and  diffusion  within  strong  magnetic  fields,  are  needed  to  perform  the 
required  flow  simulations.  The  MFD  equations  can  introduce  severe  numerical  simulation 
difficulties.  The  goal  of  this  research  is  to  develop  algorithms  for  their  solution  for  weakly 
ionized  aerodynamic  flows  for  both  internal  and  external  MFD  flows. 

The  algorithms  presented  herein  contain  numerical  procedures  for  solving  the  governing  MFD 
equations,  applying  boundary  conditions,  evaluating  temperature  and  pressyres  for  an  open 
ended  set  of  fluid  species  composing  the  flow,  and  determining  the  electrical  conductivity  of  the 
ionized  gas.  Applications  relevant  to  flows  past  flight  vehicles  and  through  scram  jet  engines  are 
presented. 


I.  Background 

There  are  several  descriptions  of  MFD  flows,  each  with  their  own  assumptions  and  set  of 
governing  equations.  The  full  set  of  MFD  equations  and  a  simplified  set,  called  the  “Low 
Magnetic  Reynolds  Number  Approximation  equations,  are  presented  below.  The  latter  set  is 
prevalent  today.  This  is  because  of  its  simplicity  and  belief  that  it  sufficiently  describes 
aerodynamic  flows  within  electromagnetic  fields.  However,  an  example  is  also  presented 
showing  that  this  may  not  be  true.  The  research  of  this  contract  was  aimed  toward  the 
development  of  algorithms  for  the  full  set  of  MFD  equations. 

The  Equations  of  Magneto-Fluid  Dynamics 

.  dU  dF  dG  dH  A  . .  TJ  r  v 

1)  1  he  Navier-Stokes  equations - h - 1- - 1- - =  0,  with  U  =  \p,pu, pv,pw,e\  , 

dt  dx  dy  dz  1  J 
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density  p ,  velocities  u,  v  and  w ,  total  energy  per  unit  volume  e . 


dE  1 

2)  Maxwell’s  equations  -  The  Ampere-Maxwell  equation  —  =  — 

dt  £e 

dB 


Faraday’s  equation  —  =  -V  x  E 


with  constraints 


VE  =  —  pc  and  V-fi  =  0, 


with  electric  Field  E ,  magnetic  field  B ,  current  density  J ,  inductive  capacity  ee , 


magnetic  permeability  pe  =  An x  10 


,-7 


kg-m 


and  charge  density  pc . 


(coulomb)2 

3)  Generalized  Ohm’s  law  J  =  cre(E  +  ux  B) ,  with  electrical  conductivity  cr , . 


1  * 


4)  MFD  assumptions:  charge  neutral  plasma,  pc  =  0,  and  time  invariant  E ,  or  J  =  — Vx/T 

Me 


Because  of  the  MFD  assumptions  in  (4),  the  above  set  of  equations  are  not  technically  the  “full” 
set  of  MFD  equations  and  we  will  respect  this  condition  by  placing  an  asterisk  on  the  word 
“full*”  used  herein.  These  assumptions  essentially  remove  the  Ampere-Maxwell  equation  from 
the  full  set  of  MFD  equations. 

The  Low  Magnetic  Reynolds  Number  Approximation 

An  ionized  flow  within  an  imposed  magnetic  field  can  self  induce,  thus  changing  the 
magnitude  of  the  total  magnetic  field.  The  relative  magnitude  of  the  induced  component  depends 
upon  the  Magnetic  Reynolds  number,  defined  by  Rm  =  u0l0cri.pe ,  where  u0  and  /0  are  reference 

flow  speed.  The  magnetic  diffusion  coefficient  is  given  by  ve  = — J — .  For  most  aerodynamic 

aeMe 

flows  the  gas  conductivity  is  very  small  and,  consequently,  ve  is  very  large  and  the  magnetic 

Reynolds  number  is  less  than  one.  In  such  cases,  any  self  induced  magnetic  field  supposedly 
rapidly  diffuses  away,  leaving  only  the  imposed  magnetic  field.  This  leads  to  a  great 
simplification  in  calculating  the  electro-magnetic  effects  upon  an  ionized  flow.  This  approach, 
the  low  magnetic  Reynolds  approximation  approach  ,  removes  the  Faraday  equation  from  the 
MFD  equation  set  and  only  needs  to  add  a  source  term  to  the  usual  flow  equations  to  include  the 
electro-magnetic  field  effects.  This  source  term  consists  of  the  Lorentz  force,  L /  =Jx  B ,  acting 
on  the  momentum  of  the  flow,  Joule  heating,  caused  by  the  flow  of  electric  current  through  the 
fluid,  plus  the  magnetic  force  work  terms  affecting  the  energy  of  the  plasma.  The  equations 
become 

I  fill  y 

—  +  —  +  —  +  —  -5,  with  S  =  \o,(JxB)x,(jxB)y,(JxB),±J  J  +  (JxB)  uy 
dt  ox  dy  dz  L  '  J 


Why  the  Need  for  the  Full*  Set  of  MFD  Equations? 
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C.  Park,  D.W.  Bogdanoff  and  U.B.  Mehta 
presented  a  1  -D  analysis  of  the  performance 
of  a  scramjet  propulsion  system 
incorporating  the  MFD  energy  bypass 
concept1.  The  MFD  accelerator  section  was 
a  square  converging  duct,  2.846m  long,  of 
height/width  0.933m  at  the  entrance  and 
0.730m  at  the  exit,  part  of  which  is  sketched 
in  Fig.  1 .  It  is  located  just  down  stream  of  the 
combustor  section.  Combustion  was  not 
simulated  numerically  in  this  study. 

The  magnetic  field  across  the  channel  was 


2  - 


1  - 


combustor 


MFD  accelerator 


j 


Figure  1 

simulation 


Sketch  of  MFD 


x  2 

accelerator 


50  =  11.28T  and  the  transverse  voltage  gradient 


varied  from  30,990  V/m  at  the  entrance  to  31,470  V/m  at  the  exit.  At  the  entrance  the  pressure 
was  1 .25 1  x  1 06  N/m\  temperature  3583°  K,  and  the  Mach  number  equaled  1.147.  The  electrical 
conductivity  was  crt,  =35.87 mho/m.  This  value  was  used  in  the  Park,  Bogdanoff  and  Mehta 


study.  The  interaction  parameter  Q-o’l,B^l0/pxux=20  and  /?m=0.17,  based  on  accelerator 

channel  length.  The  velocity  vectors  within  the  combustor  and  accelerator  sections  are  shown  in 
Figs.2  and  3.  They  both  show  acceleration  in  accord  with  the  load  factor  set  at  2.04.  Flowever,  it 
can  also  be  seen  that  the  velocities  are  larger  in  the  combustor  section  for  the  low  magnetic  Re 
number  approach. 


Figure  2  Velocity  vectors,  low  magnetic  Re  approach  Figure  3  Velocity  vectors,  full  MFD  approach 


Figure  4  Thrust  and  velocity,  low  magnetic  Re  approach  Figure  5  Thrust  and  velocity,  full*  MFD  approximation 

A  disturbance  that  has  propagated  upstream  against  the  Mach  1.147  flow  can  be  observed  in 
Figs.3  and  5.  This  can  only  be  a  shock  wave,  which  could  cause  serious  engine  unstart  and  was 
missed  by  the  low  magnetic  Reynolds  number  approach. 
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II.  Technical  Progress 


Numerical  Procedures 

The  following  six  sections  represent  numerical  procedures  that  have  been  devised  or  further 
developed  during  the  course  of  the  present  contact  period.  We  start  by  presenting  again,  though 
in  greater  detail,  the  Full*  Set  of  MFD  equations  and  the  Low  Magnetic  Reynolds  Number 
equations.  An  additional  description  that  avoids  numerical  precision  difficulties  in  flow 
simulations  within  strong  magnetic  and  electric  fields  is  also  described. 

(1)  The  Equations  of  Magneto-Fluid  Dynamics 

(i)  The  Full*  Equations  of  Magneto-Fluid-Dynamics 

The  unsteady  equations  of  compressible  viscous  flow  within  an  imposed  magnetic  field  become 


dU  dF  dG  dH  n 
dt  dx  dy  dz 

The  state  flux  vector  is  given  by 

U  =  [p,pu,  pv,  pw,  e,Bx,  By ,  5.  ] r 

with  density  p,  velocities  w,  v  and  w,  total  energy  per  unit  volume,  including  magnetic  field 
energy,  e’=e  +  -^B2,  and  Bx,  By  and  Bz  are  the  components  of  the  magnetic  field. 

e  =  p{e  +  j(u2  +  v2  +  w2)),  e  represents  the  internal  energy  and  B2  =  B2  +  B2  +  B2 .  The  flux 
vector  F  becomes 


F  = 


pu 


pu2  +  p  +rxx-ji:BxBx 


pvu  +  r 


xy 


-±B  B 

pe  x  y 


pwu  -f  t  -^-BB, 

XZ  Pt  X  4. 


(e  +  p"  +  T„  )u  +  TvV  +  zx:w~k^ 

+i :(-(B  •  m  +  M  +  PxyBy  +  Px:B: ) 

Byu-Bxv  +  /3xy 
B.u  -Bxw  +  (3x: 


with  p  =p  +  ^-B2  and  magnetic  stress  given  by  =-ve 


dB,  dB, 
ydx,  dxu 


The  magnetic  field 


components  shown  above  represent  the  total  of  the  imposed  and  induced  fields.  The  viscous 
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stress  tensor  is  given  by  ry  =  -p  ^-+ — —  where  SfJ  is  the  Kronecker  delta.  The 


other  flux  vectors  G  and  H  are  similar.  The  eigenvalues  are  for  the  flux  vector  F  are 


(ii)  The  Governing  Equations  Under  the  Low  Magnetic  Reynolds  Number  Approximation 

An  ionized  flow  within  an  imposed  magnetic  field  can  self  induce,  thus  changing  the 
magnitude  of  the  total  magnetic  field.  The  relative  magnitude  of  the  induced  component  depends 
upon  the  Magnetic  Reynolds  number,  defined  by  Rm  -  u0l0crL.pe ,  where  u0  and  l0  are  reference 

flow  speed  and  length,  crt,  is  the  gas  electrical  conductivity  and  ,  pc  =  4^xl0-7  is  the  magnetic 
permeability.  The  magnetic  diffusion  coefficient  is  given  by  ve  - — ! — .  For  most  aerodynamic 


flows  the  gas  conductivity  is  very  small  and,  consequently,  ve  is  very  large  and  the  magnetic 

Reynolds  number  is  less  than  one.  In  such  cases,  any  self  induced  magnetic  field  supposedly 
rapidly  diffuses  away,  leaving  only  the  imposed  magnetic  field.  This  leads  to  a  great 
simplification  in  calculating  the  electro-magnetic  effects  upon  an  ionized  flow.  This  approach, 
The  Low  Magnetic  Reynolds  Approximation  approach,  only  needs  to  add  a  source  term  to  the 
usual  flow  equations  to  include  the  electro-magnetic  field  effects. 

Ionized  flow  in  the  presence  of  a  magnetic  or  electric  field  generates  a  current,  according  to 
Ohm's  law,  j  =  cri,(E  +  uxB),  where  j  is  the  current  density,  E  is  the  electric  field  potential,  w 

is  the  flow  velocity  and  B  is  the  magnetic  field.  The  electric  current  itself  interacts  with  the 
magnetic  field  to  create  a  Lorentz  force,  Z,/  =  jx  B,  that  acts  on  the  flow  in  addition  to  pressure, 
p,  and  viscous  stress. 

In  addition  to  the  Lorentz  force  added  to  the  momentum  equations,  Joule  heating,  caused  by 
the  flow  of  electric  current  through  the  fluid,  plus  magnetic  force  work  terms  need  to  be  added  to 
the  energy  equation.  The  equations  become 


dU  dF  dG  dH 
dt  dx  dy  dz 


The  state,  flux  vectors  and  the  source  vector  are  given  by  U  =  [p,  pu ,  pv,  pw,e\ 
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> 

r  „  \ 

pu 

0 

pu'+p  +  T^ 

(7x5), 

pVU  +  Txy 

,etc,  and  S  — 

(7x5), 

pwu  +  rX2 

(7x5), 

,  X  7  dT 

(e  +  p  +  Txx)u  +  TxyV  +  Tx:w  -k  — 
V  OX  J 

v(7x5)«+j-7-7y 

The  eigenvalues  for  the  flux  vector  F  ,  under  the  low  Magnetic  Reynolds  approximation,  are 

\ 3 4  =  u  and  X15=u±c 

This  equation  set  is  supposedly  sufficient  to  describe  ionized  flow  within  an  electro-magnetic 
field,  as  long  as  the  fields  are  specified.  It  is  not  much  more  difficult  to  solve  than  the  underlying 
Navier-Stokes  flow  equations  themselves.  However,  if  the  magnetic  field  varies  in  time  by  self 
induction  and  the  induced  magnetic  components  are  relatively  significant  in  magnitude  then  the 
equations  for  magnetic  induction  also  need  to  be  solved.  This  larger  set,  shown  earlier,  is  much 
more  difficult  to  solve  and  should  not  be  attempted  if  it  can  be  avoided.  However,  as  shown  in 
the  Background  section  above,  this  simpler  set  of  equations  may  not  suffice  even  for 
aerodynamic  flows. 

(in)  An  Alternative  Formulation  of  the  Equations  of  Magneto-fluid-Dvnamics  - 

The  Reduced  MFD  Equations 

An  alternate  formulation  of  the  governing  equations  of  magneto-fluid-dynamics  has  been 
devised.  This  formulation  is  mathematically  equivalent  to  the  original  set  of  equations  governing 
magneto-fluid-dynamics,  given  above  in  Sec.(i),  including  magnetic  self  induction,  and  retains 
the  conservation  law  form  of  the  equations  and  their  eigenvalues.  This  new  set  has  advantages 
over  the  original  set  when  solved  numerically  for  flows  within  strong  imposed  magnetic  fields. 

This  formulation  treats  the  imposed  and  induced  fields  separately.  Though  ever  present,  the 
imposed  field  remains  in  the  background  as  the  finite  volume  difference  equations  focus  on  the 
induced  field.  The  imposed  field  can  not  be  eliminated  entirely  from  the  difference  equations 
because  of  non-linearity,  but  no  squared  terms  of  the  imposed  magnetic  field  appear.  This  is 
important  for  flows  within  strong  imposed  magnetic  fields,  because  the  magnetic  pressure, 
proportional  to  the  square  of  the  imposed  magnetic  field  strength,  can  be  several  orders  of 
magnitude  larger  than  the  aerodynamic  pressure  or  the  induced  magnetic  field  pressure. 
Numerical  errors  in  the  very  large  magnetic  stress  difference  terms  of  the  imposed  field  could  be 
of  significance  when  combined  with  the  relatively  smaller  fluid  stress  terms. 

Before  the  introduction  of  this  alternative  formulation,  there  were  two  choices  for  including 
the  effects  of  a  magnetic  field  upon  an  ionized  flow:  (1)  the  complete  equations  of  magneto¬ 
fluid-dynamics,  including  magnetic  induction  and  diffusion  and  (2)  the  inclusion  of  the  “j  cross 
B”  force  and  Joule  heating  effects  in  the  Navier-Stokes  equations  as  additional  source  terms.  The 
first  choice  is  a  set  of  eight  equations  consisting  of  the  Navier-Stokes  equations,  with  added 
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magnetic  stress  tensor,  plus  the  Maxwell  inductions  equations  and  is  presented  in  Sec.(i)  above. 
The  second  choice,  a  set  of  five  equations,  called  the  “Low  Magnetic  Reynolds  Number 
Approximation”,  assumes  that  the  induced  magnetic  field  is  negligible.  It  is  far  more  efficient 
and  has  far  fewer  numerical  difficulties  associated  with  inclusion  of  magnetic  effects  than  the 
first  choice.  There  is,  however,  some  uncertainty  in  when  the  Low  Magnetic  Reynolds  number  is 
valid,  even  for  aerodynamic  flows  of  current  interest.  The  now  available  third  choice  presented 
below  is  a  consistent  alternative  to  the  first  choice. 


The  total  magnetic  field  consists  of  the  imposed  magnetic  field  Bo  and  the  induced  magnetic 

field  Bi .  Bi  =  Bo  +  Bi ,  where  the  subscripts  t,  0  and  /  now  and  below  indicate  total,  imposed  and 
induced  magnetic  components.  For  cases  for  which  the  induced  field  is  much  less  than  the 
imposed  field,  but  not  negligibly  small,  we  can  benefit  by  rewriting  the  Lorentz  force  as 

Lf  =  —  (V x  B,)x  Bi ,  because  the  imposed  magnetic  field  is  generated  by  currents  external  to 

the  flow  field,  for  which  VxBo-0.  The  approach  taken  here  is  similar  to  the  simplification  in 
electromagnetic  scattering  where  only  the  disturbed  field  is  calculated,  with  nothing  lost  by  the 
separation  of  the  two  fields. 

We  can  also  write  the  Lorentz  force  as  Lf  =  —  (V x  Bi)x  Bi  -  — (V x  Bo)x  Bo,  and  through 

Me  Me 

some  algebraic  manipulation,  the  Lorentz  force  can  be  brought  into  the  flux  derivative  terms  of 
the  momentum  equation,  in  conservation  law  form  as  before.  The  state  vector  becomes 

U  =  [p,pu,pv,p\v,e',B,x,Biy,Bi .]  and  the  new  flux  vector  F  becomes 


F  = 


pu 

pu1  +  P  +  -j;BixBlxVV-j-{BixBox} 

PVU  +  Txy  —kBiyB'x  -k\B'xB°y\ 

pwu  +  Tx:  -jrB-.B,  -j;{B,xBo:) 

....  .  dT 

(e  +p  +Txx)u  +  T xyv+Tx.w-k  — 

ox 


+  — 


S)B,,  +  P„B,,  +  PVB,,  +  flA ) 


B<yu-B,xv  +  pxy 
Bi.u  -  Bixw  +  /?xz 


with  p’  =  p  +  -±-Bi2  +yB,  Bi  and  the  magnetic  stress  given  by  (3tl  =  —ve 

By  replacing  the  magnetic  pressure  B r  by  the  smaller  y '—B,2+-^B,  Bi,  the  magnetic  and 

static  pressures  are  closer  in  magnitude  for  strong  imposed  magnetic  fields.  Favorable  reductions 
also  take  place  in  the  induction  equations  because  the  magnetic  diffusion  terms  y9(/  are  just 


SB,,  SBi, 


dx.  dx. 


7 


components  of  the  curl  of  the  induced  magnetic  field  times  ve .  Again  the  imposed  field, 
produced  by  currents  outside  the  flow  field,  is  curl  free.  Hence,  vtVxfi  =  v(Vxfi.  Here  also 

the  production  and  diffusion  terms  are  more  equally  balanced.  Finally,  the  terms  in  the  curly 
brackets  in  the  momentum  equations  above  vanish  if  the  imposed  field  is  constant  in  space 
because  of  the  divergence  free  nature  of  both  the  imposed  and  induced  fields. 

One  may  assume  that  the  structural  changes  to  the  equations,  just  presented,  from  the 
separation  of  the  induced  and  imposed  fields,  would  have  profound  changes  to  the  original 
eigenvalue-eigenvector  structure  of  the  equations.  Fortunately,  the  eigenvalues  remain  the  same 
as  shown  below. 


B, 


4,6="’  4m  =m±TJ,/12.5  =u± 


f 

1  . . 

2  — 'A 

1 

«*+£‘  +  J 

-4c2 

2 

P  \ 

p  J 

P 

and 


^7.8  =U± 


: 2  ,  B, 


2  ,  B, 


C  +—  -J  c  +  — 


p  J 


-4c2  — 


The  eigenvectors  are  changed,  however,  but  the  original  set  can  still  be  used  in  the  solution 
procedure  as  is  to  solve  the  alternative  RMFD  (Reduced  Magneto-Fluid  Dynamics)  equations 
just  presented,  because  of  the  conservation  form  of  the  flux  vector  splitting  used. 


The  term  Reduced  is  used  here  to  reflect  the  notion  that  the  magnitude  of  the  magnetic  terms 
is  reduced  by  removing  as  often  as  possible  the  imposed  magnetic  field  from  them,  although  the 
number  of  terms  is  actually  increased.  The  RMFD  equations  are  mathematically  and  physically 
equivalent  to  the  original  MFD  equations. 


(2)  Boundary  Conditions 

til  At  Solid  Walls 

The  usual  “no  slip”  type  boundary  conditions  are  used  for  the  Navier-Stokes  terms  of  the 
governing  equations.  The  boundary  conditions  for  the  magnetic  terms  appearing  in  the  fluid  flow 
equations,  i.e.,  the  fluid  momentum  equations  containing  the  Lorentz  force  terms  and  the  energy 
equation  containing  the  Joule  heating  and  magnetic  force  work  terms,  are  as  follows. 


a)  The  imposed  magnetic  field  is  specified  everywhere. 

b)  The  normal  component  of  the  induced  magnetic  field  is  fixed  to  zero  (see  note  below). 

c)  The  tangential  components  of  the  induced  magnetic  field  are  chosen  so  that  V  x  B,  =  0 . 

d) 

f  r,  \ 

,  £is  specified  at  the  wall. 


SB 


For  the  induction  equation,  —  =  -V  x  E  =  -V  x 


-  R  — 

-uxB-\ - Vx  B 


Note  that  for  perfectly  conducting  walls  (cr,  =  oo)  and  the  “no  slip”  boundary  condition  implying 
that  u  =  0,  the  tangential  components  of  E  must  vanish,  otherwise  the  current  J  a  a  E  would 
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be  infinite.  This  implies  that  the  normal  component  of  — must  vanish,  proving  assertion  (b) 

dt 

above. 

(ii)  At  Symmetry  Planes.  Entrance.  Exit  and  Far  Field  Boundaries 

At  planes  of  flow  symmetry  reflection  boundary  conditions  are  applied.  At  flow  through 
boundaries  the  eigenvalues  of  the  equations  should  be  used  to  determine  the  dependence  of  the 
boundary  on  the  interior  of  the  flow.  The  equations  of  MFD  have  characteristic  speeds  that  differ 
greatly  from  the  wave  speeds  of  the  conventional  equations  for  compressible  flow.  Boundaries 
that  appear  to  be  supersonic  may  be  actually  subsonic  boundaries  because  of  fast  magnetic 
waves.  This  is  a  significant  concern  for  using  the  “Low  Magnetic  Reynolds  Number 
Approximation”  approach,  which  is  blind  to  this  boundary  condition  dependence  issue. 

(3)  Electrical  Conductivity  and  Ohm’s  law 

Most  previous  numerical  studies  for  solving  the  MFD  equations  used  a  scalar  electrical 
conductivity.  However,  a  scalar  electrical  conductivity  can  not  simulate  the  effects  of  Hall 
currents  and  ion  slip.  A  tensor  form  for  the  electrical  conductivity  of  the  gas  is  required  and  is 
described  below. 

(i)  Scalar  Conductivity 

The  current  density  J can  be  defined  by  Ohm’s  law,  using  a  scalar  electrical  conductivity  cr, ,  as 
follows 

J  =  cre(E  +  mx  B) 

The  conductivity  depends  upon  the  number  and  mobility  of  the  charged  particles,  both  electrons 
and  ions,  present.  For  air  it  is  usually  very  small  and  consequently  the  magnetic  and  electric 

fields  required  for  aerodynamic  interaction  need  to  be  very  large.  If  we  define  E'  =  E  +  w  x  B  , 
then  current  density  can  be  written  as  J  =  acE' .  Thus  J  is  in  the  same  direction  as  E' .  This  is 
not  strictly  true  because  particle  collisions  also  can  cause  the  charged  particle  to  drift  off  this 
course,  the  E'xB  drift,  within  a  plane  normal  to  B ,  which  causes  the  Hall  current  and  ion  slip 
effects  and  requires  the  electrical  conductivity  to  be  represented  as  a  tensor. 

(ii)  Tensor  Conductivity  (following  along  the  discussion  in  Mitchner  and  Kruger,  Partially 
Ionized  Gases') 

Consider  a  local  coordinate  system  (x,y,z),  where  z  is  in  the  direction  of  B,  jp  is  in  the 
direction  of  E'  in  the  plane  1  to  B  and  x  is  J.  to  both  y  and  z .  Then  B-Bis,  and 
generalized  currents  in  the  x-y  are  denoted  as  J-  =  <JLE-y  and  Ji=crHE'i,,  and  in  the 
£  direction  by  J.  =anEi.  Consider  a  second  coordinate  system  (. x\y\z '),  which  represents  a 
coordinate  rotation  about  the  £  axis  through  a  angle  6 .  Then 

Jx,  =  cos  0J-X  +  sin  OJ- 
J„.  =-sin6L/.  +cos#J- 

y  x  y 

Or 
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Jx.  =  -crH  cos  9E -  +  crx  sin  9E-  =  -<JHEy.  +  <r±Ex. 

Jy,  =  <th  sin  #£.  +  <tx  cosOE-  =  crHEx.  +  cr1£l,. 

Or,  by  defining  the  electrical  conductivity  tensor  E,  the  current  density  J  ,  via  Ohm's  law,  can 
be  written  as 


Jx- 

Jy 

Js, 


>1 

0  ' 

~K~ 

0 

E'y 

0 

0 

a« 

E2, 

E 


The  first  coordinate  system,  ( x,y,z )  is  aligned  with  the  local  magnetic  and  electric  field 
directions.  The  second  ( x’,y',z ')  generalizes  the  direction  within  the  y-z  plane.  We  now  need 
to  relate  these  to  the  computational  coordinate  system  (x,y,z) .  We  define  again  the  (x,y,z) 

-  _  B  E, 

coordinate  directions  as  follows  z  -b  =  -, — - .  E„  =  (E'b)b  ,  E.=E'-E„,  y-e=rA  and 

|5|  "  y  "  |£.| 

x  =  yxz  =  d  .  Then,  using  the  orthogonal  unit  vectors  b  ,  e  and  d defined  in  the  ( x,y,z ) 


coordinate  system,  b  =  (bx,by,b2),  e  =  (ex,ey,e.)  and 


-  yux ,  w  y9  u,  ) 


Or 


Similarly, 


X 

dx  dy 

d. 

X 

X 

dx  ex 

K' 

X 

y 

= 

e, 

y 

and 

y 

= 

d>  ey 

by 

y 

z 

0 

.  Sr- 

* 

z 

z 

A  e2 

z 

T 

T~ ' 

x'~ 

cos  0  sin# 

0] 

X 

X 

cos  6 

-sin#  0 

y' 

= 

-sin#  cos# 

0 

y 

and 

y 

= 

sin# 

cos#  0 

z'_ 

_ 

0 

0 

1 

z 

z 

0 

0 

lj 

- - '  ' - V— 

5  S"1 


Starting  with  the  current  density  defined  above  J  =  E£\  written  in  the  (x\y\z')  coordinate 

system,  we  can  express  J  in  the  computational  coordinate  system  (x,y,z),  using  the 
transformations  S  and  T ,  as  follows. 


~JX~ 

T~'S~' 

Jy 

=  r'S-'l.STT-lS-' 

E'y 

Jr. 

K 
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Or 


v/ 

Ex+(uxB\ 

Jy 

=  r'S~'l.ST 

K 

=  R 

Ey+(ux  B)y 

J. 

R-' 

E. 

E.  +(ux  B). 

~z~ 

1 

+ 

/> — s 

Si 

X 

&01 

x 

_ 1 

K 

- 

Ey+(ux  B)y 

=  T 

V - v - * 

Jy 

E\ 

E.  +  (u  x  B). 

R 

y,_ 

,  where  I  = 


crl+a2H 


1 


°T  °h  0 


_2  .  2 

0  o 


and  /?  is  the  resistivity  matrix.  Using  the  MFD  assumption,  J  =  — V  x  B ,  we  can  express  the 

Me 

electric  field  as 

E  =  -uxB  + —  VxZ? 


The  equation  for  magnetic  induction  is  then 


^,-Vx£  =  -Vx 
dt 


f  ^  /?  —  __  ^ 

-uxB-\ - Vx5 


v 


/ 

For  the  scalar  electrical  conductivity  R  =  —  ,  otherwise  R  =  T  .  The  conductivities 

°e 

used  within  the  tensor  L'1  are  defined  by 

1  +  5 


with 


(\  +  sf+p2 


<3e >  <JH  — 


Pe 


2  ,  ril  e 


(i +sy+pt 


<?e  and  a „  =  cre  (1) 


s  =  the  ion  slip  factor  for  a  weakly  ionized  gas. 


CO 


Pe  =-±,  the  Hall  Parameter  for  electrons. 


co, 


Pt  =-3-,  the  Hall  Parameter  for  single  charged  heavy  ions. 


v, 


The  cyclotron  (or  Larmor)  frequency  for  the  electrons,  in  terms  of  the  magnetic  field  strength, 
electron  charge  and  mass,  is  given  by  coe  =  the  corresponding  frequency  for  the  ions  is 


m. 


co.  = 


kl* 


m. 


,  the  average  momentum  transfer  collision  frequency  for  electron  is  ve  and  that  for 


ions  is  v( .  For  the  present  flow  simulations  /?,  =  0.5  and  pt  =  0.05 
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(4)  Equilibrium  Chemistry  Model  and  Evaluation  of  cr. 

Electrons  and  ions  are  created  by  collisions  between  the  atoms  and  molecules  composing  the  gas, 
in  our  case  air  and  any  element  introduced  into  the  flow  via  seeding.  The  collisions  become  more 
violent  with  increasing  temperature.  Thus,  ionization  depends  upon  the  chemical  mixture  of  the 
gas  and  temperature.  The  air  chemistry  model,  with  cesium  seeding,  consisted  of  1 3  species  and 
1 1  reactions.  The  species  consisted  of 

N2,  02,  NO,  N,  O,  Cs,  e~ ,  N2, 0*2  ,  NO+ ,  N+ ,  0+  and  Cs+ 

The  reactions  1 1  considered  are 

N2+M  ++2N  +  M 
02  +  M  ++20  +  M 
NO  +  M  ++N  +  O  +  M 
N2+0  ++NO+N 
NO  +  O  ++N  +  02 
N  +  O  NO*  +  e~ 

For  the  cases  to  be  presented  cesium  seeding  was  not  used.  The  gas  will  be  assumed  to  be  in 
chemical  equilibrium.  Equilibrium  coefficients,  functions  of  temperature  and  particle  number 
densities  which  determine  specie  concentrations,  were  taken  from  Park’s,  Nonequilibrium 
Hypersonic  Aerothermodynamics3 .  The  local  mass  fractions  were  then  used  as  inputs  to  Park’s 
program  “tmpt”4  for  calculating  transport  coefficients,  including  cr, . 

Equilibrium  Calculation 

The  set  of  species  considered  above  is  an  open  set  in  that  a  new  species  can  be  added  at  will. 
This  is  different  than  the  usual  fixed  set  of  species  with  temperature  and  pressure  determined 
from  a  set  of  curve  fits  for  the  given  fixed  set  of  species.  However,  temperature  and  pressure 
must  be  calculated  instead  from  chemical  and  thermal  equilibrium  relations,  which  is  more  time 
consuming.  Considerable  effort  has  been  made  to  compute  equilibrium  conditions  efficiently. 
The  flow  solver  for  the  governing  MFD  equations  solves  the  conservation  laws  for  mass, 
momentum  and  energy.  This  determines  the  fluid  density  and  internal  energy  at  each  point  of  the 
computational  mesh.  Temperature  and  pressure  are  then  determined  from  the  fluid  density  and 
internal  energy  at  each  point  by  solving  a  set  of  equilibrium  relations.  The  species  are  ordered  by 

{5,. }  =  { N2,02,  NO,  N0*,N2,02,N,  0,e-,Cs,Cs+,N+,0+ } 

The  procedure  for  calculating  pressure  p  and  temperature  T ,  given  density  p  and  internal 
energy  e ,  at  each  grid  point  is  as  follows. 

1)  The  set  of  concentrations  is  initialized,  c*  =c,°  = /?“  /  p,  i  =  \,-,ns,  where  ns  is  the 
number  of  species,  p°  and  c(°  are  the  initial  density  and  concentration  of  species  i  and 
the  superscript  k  indicates  the  newest  value.  These  initial  values  correspond  to  those  for 


Cs  *+■  e  ++  Cs  +  e  +  e 
N  +  e~  <->  N+  +  e~  +e~ 
0  +  e~  <->  O*  +e~  +  e~ 
N  +  N  ++N2+e' 

O  +  O  ++02+e~ 
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a  given  reference  temperature.  For  example,  for  a  reference  temperature  300°  K  with  no 
cesium  seeding , 

c®  =  0.763,  c\  =  0.237,  c“>2  =  0 

The  reference  temperature  of  300°  K  was  chosen  in  the  present  calculations.  A  better  choice 
would  be  the  values  for  c,  from  the  previous  time  step,  but  this  would  require  significant 
storage  memory,  which  is  avoided  at  present. 

2)  The  temperature  T  is  initialized  by  the  previous  temperature  from  the  last  time  step  at  the 
given  grid  point  and  is  given  by  jk-j° 

3)  The  set  of  concentrations  jc*  j  consistent  with  temperature  Tk  is  found  by  solving  a  set  of 

equilibrium  relations.  The  number  of  kg-moles  of  species  i  per  unit  volume  is 
X,  =  p,  / M: ,  where  Mt  is  the  molecular  weight  in  kg  of  species  /.  The  chemical 
reactions  listed  above  are  of  two  types 

I.  Dissociation:  molecule  AB  breaks  apart  to  form  species  A  and  B,  AB  <-»  A  +  B 
The  chemical  reaction  relation  is  k'f xAB  <-»  k'h xAXn ,  where  k'f  and  k'h  are  the  forward 
and  backward  reaction  rates  for  chemical  reaction  AB  <-»  A  +  B .  In  equilibrium  this 

k' 

relation  becomes  klfxAB=k,bzAZB  or  kL,Xab=ZaZb>  where  KL,  =  T7  The  forward 

and  backward  reaction  rates  were  obtained  from  tabulated  values  in  Park3  .  They 
depend  upon  both  particle  number  density  and  temperature. 

II.  Exchange  reaction:  species  A  and  B  form  species  C  and  D,  A  +  B  <-»  C  +  D . 
Similarly,  the  equilibrium  relation  for  this  reaction  is  K'^XaXb  =  XcXd  ■ 

Qnly  quadratic  equations  need  to  be  solved  for  each  reaction  equation  to  bring  the 
concentrations  into  equilibrium.  Consider  that  an  additional  x  kg-moles/m3  need  to  react 
to  bring  about  equilibrium.  For  a  type  I  chemical  reaction  the  equation  is 
k'<,(Xab-x)  =  (Xa+x)(X b+x)>  and  for  type  II  K'^x/I-x)(xh-x)  =  (xc+x)(xd+x). 
Positive  x  indicates  that  the  reaction  is  “forward”.  Checks  need  to  be  made  so  that  the 
species  concentrations  remain  non-zero.  For  type  I  reactions  .v  needs  to  bounded  by 
-mi x\(xa,Xb)-x-Xab  If  A-B  then  also  -\xA-x-  F°r  type  II  reactions 

-min(^c,  jD)<x:<min(^,^B)  .  But  if  A-B ,  then  also  x<\xA  and  if  C  =  D,  then 
also  Xc  -  x  ■  If  x  *'es  outside  bound  limits  it  is  reset  to  the  nearest  limit. 

4)  The  checking  of  the  kg-mole/m3  change  x  for  each  reaction  in  the  above  step  is  still  not 
sufficient  to  avoid  negative  concentrations.  Each  reaction  was  checked  independently  of 
the  others  and  when  combined  may  still  lead  to  negative  concentrations.  A  global 
checking  needs  to  be  made.  For  example  species  N2  appears  in  the  first  and  fourth 

reactions  in  column  one  of  the  listing  above.  If  x(1)  and  x(41  indicate  changes  calculated 
for  each  reaction,  then  the  kg-mole/m3  change  to  N2  needs  to  be  bounded  by 

x(l)  +x(4)  =  xw  <  Xn2  •  To  insure  that  all  the  concentrations  remain  non-negative  we  first 
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set  r}  =  1,  j  =  l,-,nr,  where  nr  is  the  number  of  reactions  (i.e.,  nr  =  11  for  the  above 

set  of  chemical  reactions)  and  denote  x(J),j  =  1  ,  -,nr  as  the  set  of  kg-mole/m3  changes 
for  the  set  of  chemical  reactions.  Some  of  the  reactions  need  to  be  limited  to  avoid 
negative  concentrations  from  occurring.  Each  species  needs  to  be  checked  in  turn.  For 
example,  for  the  first  species,  if  x(l)+x<4)  =  xNi  >  xNl ,  then  the  reactions  involving  N2 

x(')+x<4)  *0)  +  (4) 

are  limited,  as  follows.  r{  <-  min(rp - - — )  and  r4  min(r4,: - ) .  Similarly, 

Xn2  X  n2 

for  the  second  species  02,  which  appears  in  reactions  (2)  and  (5),  if 
xr.)  _JC(5>  _  >  (notice  the  minus  sign  appearing  before  x<5),  which  results  from  its 

x<2)  -  v(5) 

presence  on  the  “backward”  side  of  the  reaction  equation),  then  r2  <—  min(A-,,- - - — ) 

Xo2 


-x(5) 

and  r5<-min(r5, - ),  etc.  After  checking  each  species  in  turn  the  set  of  kg- 

Xo2 

mole/m3  changes  is  reset  x(j)  <-rjXU),  j  ~\,---,nr  and  each  species  is  adjusted  toward 
equilibrium,  for  example,  XNi  <-  XNi-xw -xw  and  Xo2  Xo2  ~ *(2)  +  *(5) »  etc-  The 


new  concentrations  become  c*+l  <-  aM,  —  +  (1  -a)ck ,  i  =  l,---,ns,  where  a  is  a 

P 

relaxation  parameter,  currently  taken  to  be  V2. 

5)  The  new  concentrations  will  distribute  the  internal  energy  differently  and  the 

ns  R 

thermodynamic  properties  of  the  gas  need  to  be  recalculated,  as  follows.  R  =  ^c, 


l 

M 


where  R  and  R  are  the  specific  and  universal  gas  constants  and  the  superscript  k  on 


the  concentrations  has  been  surpressed;  cv  =^jcl 


3  Rg 

2  A 4,  +C“ 

\  irans 


nds 

1 


M. 


■  +  c, 


vib- 


\  rot 


where 


cv  is  the  specific  heat  at  constant  volume,  nds  is  the  number  of  diatomic  species  (  the 
first  nds  species  in  the  ordered  list).  The  specific  heat  of  the  gas  contains  contributions 
from  translational,  rotational,  vibrational  and  electronic  state  internal  energy  and  cel  and 

cvih  are  functions  of  temperature  Tk .  The  total  internal  energy,  with  concentration  set 


ns 

{c,}at  temperature  Tk ,  equals  ek  =cvTk  +  '^jclhfl° ,  where  hf°  is  the  heat  of  formation 

1 

of  species  i  relative  to  the  reference  temperature  300°  K .  If  the  set  of  concentrations  and 
temperature  Tk  are  such  that  ek  -  e ,  then  we  have  consistency  and  we  may  proceed  to 
the  next  step  (6).  If  they  do  not  match  beyond  a  reasonable  tolerance,  we  then  need  to 
choose  a  different  value  for  Tk ,  return  to  step  (3)  and  iterate  until  convergence,  ek  -»e. 
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choice  for  the  new  estimate  of 


Tk  <—Tk  +  sign 


mm 


g*-g 

vc>'+<V 


,1000"/: 


,e  -e 


the  temperature 

ns 


IS 


,  where  chf  * 


i 


dTk 


6)  The  new  temperature  becomes  the  Tk  for  which  ek  — >  e  sufficiently.  Thus  T  =  Tk  and 
the  pressure  becomes  p  =  pRT . 


(5)  Divergence  Control  of  the  Induced  Magnetic  Field 

An  important  constraint  on  the  magnetic  field  is  that  it  be  divergence  free,  V  B  =  0 .  During 
the  course  of  the  calculation,  error  the  in  numerical  solution  of  the  magnetic  induction  equations, 
containing  both  production  terms  and  strong  diffusion  terms,  can  cause  the  magnetic  field  to  drift 
away  from  being  divergence  free.  If  uncorrected,  this  can  quickly  escalate  the  error  in  the 
simulation.  There  are  a  few  options  for  controlling  this  error:  (1)  a  convection  term  can  be 

retained  in  the  governing  equations,  proportional  to  V  •  B ,  to  carry  error  out  of  the  flow  field,  (2) 
a  poisson  equation  can  be  setup  whose  solution  can  be  used  to  drive  the  divergence  error  back 
toward  zero,  and  (3)  the  divergence  error  can  be  placed  in  a  dissipation  equation  to  attenuate  it 
toward  zero.  The  choice  of  which  procedure  to  use  is  very  case  dependent.  Procedure  ( 1 )  is  fairly 
popular,  but  causes  the  governing  equations  to  lose  the  “conservation  law”  form.  Previously,  the 
present  effort  used  procedure  (2)  for  internal  flows  with  good  success.  Procedure  (3),  perhaps 
better  suited  for  the  unbounded  external  flow,  is  believed  to  be  a  new  algorithm  improvement 
and  will  be  presented  here.  Consider  the  following  three  equations  in  “conservation  law”  form. 

8Bt  dVB  dB  dVB  SB,  dV  ■  B 

— -  = - ,  — L  = -  and  — *-  = -  (2) 

dt  dx  dt  dy  dt  dz 

If  V  B  =  0  everywhere,  the  above  equations  will  not  change  the  magnetic  field.  On  the  other 

hand,  if  V-  B  *  0,  the  equations  will  change  the  magnetic  field,  but  the  new  field  will  also  have  a 
changed  divergence  governed  by 

dVB  d2VB  d2VB  d2VB 

- — - - - 1 - - - 1 - - —  (3) 

dt  dx 2  dy 2  dz2 

This  equation  represents  the  dissipation  equation  for  the  divergence  of  the  magnetic  field.  The 
boundary  conditions  for  the  above  equations  are  that  V  5  =  0  along  all  boundaries.  The  error 
introduced  within  the  flow  field  during  a  flow  simulation  probably  has  both  positive  and  negative 

values  for  V-5,  occurring  in  equal  measure.  The  dissipation  equation  together  with  the 
boundary  conditions  should  drive  the  error  toward  zero.  However,  the  procedure  is  numerical 
and  consideration  of  how  the  derivative  terms  above  are  represented  by  finite  differences  is 
important.  Using  the  notation  for  backward,  forward  and  central  difference  operators  as  follows, 
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A- 


B  -B  B  -B 

d  x  x  n  xx 

^_1±k - ,-yj L>  =_^u> - KJ±_  and  ^o_ B 

Ax  Xij,k  x,  -  x,_,  Ax  A,*  x(+l  -x.  Ax  Xij,k  x,+1  -  x,_, 


5  -5 

jc  x 

i+l,j,k 


etc.,  where  the  “dots”  appearing  in  the  operator  notation  imply  that  the  operator  applies  to  all 
factors  to  the  right.  The  equations  above  are  differenced  as  follows.  First  the  divergence  of  the 
magnetic  field  is  approximated  by 

V-Buk  ~  —  Bx  +^BV  +^-B. 

Ax  '•JX  Ay  y,’J’k  A z 

Equations  (2)  are  then  approximated  by,  where  the  superscript  index  n  indicates  the  time  step. 


Bn+X  -B" 


A  _  B« 

AA - ■+L  =  £^V.B".,  -AiA - =  and 

A t  Ax  ’h  At  Ay 


b"+>  -b: 

~i,j,k  _  A 


At 


A z 


V  B"  . 

i,j,k 


Equation  (3),  though  not  solved,  predicts  that  the  solution  to  the  above  difference  equations  will 
result  in 


V '  B"j!k  =  V  •  5"y>i  +  A/  • 


A- A- 

Ax  Ax 


V  Bn 


ij,k 


+  D+’  D-  v  b: 


Ay  Ay 


'J,k 


A-A 

A  z  A z 


-v-b: 


i,j,k 


which  indicates  that  the  new  divergence  of  the  magnetic  field  equals  the  previous  value  modified 
by  dissipation. 

For  the  present  simulations,  Eqs.(2)  were  solved  in  a  general  curvilinear  axisymmetric 
coordinate  system.  To  preserve  axisymmetry,  central  difference  operators  were  used  to 
approximate  derivatives  in  the  z  direction. 


(6)  The  Numerical  Method 

^  .  .  dU  8F  dG  dH  A  .  ,  .  ,  , 

The  governing  equation  - + - 1-  —  + - =  0,  shown  above,  contains  both  inviscid  terms 

dt  dx  dy  dz 

and  viscous  terms.  All  terms  and  the  boundary  conditions  are  treated  implicitly.  The  inviscid 
terms,  essentially  the  Euler  equations  extended  to  the  ideal  equations  of  MFD,  use  a  modified 
Steger-Warming  flux  splitting  procedure.  Assuming,  for  the  moment  two  dimensional  flow  and 
that  the  flux  vectors  shown  below  contain  only  the  inviscid  terms,  we  can  write 


n 

F  —  AU  and  G  —  BU,  where  A  = -  and  B  = - 

dU  dU 
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We  can  turn  these  inside  out  as  follows,  assuming  also  for  the  moment  that  the  inverses  exist, 
(i.e.,  the  eigenvalues  of  A  and  B  do  not  vanish). 

U  =  A~'F  and  U  =  B~'G ,  where  A~l  =  S~'C~1A~JCaS  and  B~'  =  S-'C;'A-'CbS 

Therefore,  we  can  split  the  fluxes  F  and  G  directly  as  follows 

F±  =  S~'C~'Aa  A~'CaS  F  =  S-'CaDa  CaS  F  =  &±F 

and 

g±  =  s-'c-b'aba-'cbs  G = s-'cb'dbcbs  G  =  m±G 


where  the  diagonal  matrices  DA  and  DB ,  have  “ones”  or  “zeros”  as  elements.  For  example,  D , 

has  “ones”  where  A  A  has  non-zero  elements  and  D ,  =  I  -  DA  .  This  approach  has  the  benefit  of 

partitioning  the  conservative  flux  vectors  themselves,  which,  for  example,  unlike  the  state 
vectors,  are  often  continuous  across  discontinuities  in  the  flow.  In  terms  of  the  generic  flux 
vectors 


F*„,  =  a;  F'+a”  f" 


/+I/2.J 


i.J 


i+\.j 


MJ 


and  G"j+U2 


G"  +W  G" 


'.7+1 


i,j+l 


Higher  order  accurate  approximations  for  the  fluxes  can  be  made  by  upwind  extrapolation  or 
interpolation  of  the  flow  variables  to  the  flux  surfaces,  i.e., 


F" 

'+1/2,7 


=  &  F"  +  3"  F,"  ,  etc. 

+/+i,y  R  -/+! ,j  L 


The  method  used  herein  was  third  order  accurate  in  space  and,  because  the  flows  converged  in 
time  to  a  steady  state,  only  first  order  accurate  in  time.  The  method  is  also  TVD  (Total  Variation 
Diminishing)  to  be  discussed  below.  Also,  the  full  set  of  Navier-Stokes  viscous  terms,  plus  the 
magnetic  dispersion  terms,  veV  x  B ,  of  the  Faraday  equation,  were  included  using  second  order 
accurate  central  difference  approximations. 

TVD  (Total  Variation  Diminishing)  Method  Extension 

We  present  the  extension  to  TVD  of  the  Modified  Steger- Warming  Method.  This  extension  is 
more  than  cosmetic.  It  is  designed  to  prevent  spurious  oscillation  in  the  solution  that  can  cause 
severe  numerical  difficulties.  We  illustrate  the  procedure  by  applying  it  to  the  Euler  equations  in 
conservative  form  in  one  spatial  dimension,  as  follows.  The  Euler  equations  are 


dU  dF_ 
dt  dx 


fP  \ 

(pu  \ 

0 ,  where  U  = 

pu 

and  F  = 

pu2  +  p 

S  7 

K(e  +  p)uy 

A  generic  algorithm  for  solving  these  equations  is 
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+,  _  A/r 


Ax 


Y  i+l/2 


The  Roe  and  Modified  Steger  Warming  methods  for  evaluating  the  flux  vector  are  shown  below. 
They  are  both  first  order  methods  at  present. 


Fn 

ri+\/2 


j-y(Roe)  _  Fj  +  Fj+\ 

A+l/2  “  - 


F{M-S-W)  _  & 

/  f+1/2  ~^+i+\/2 


2 


"l"'^~/+l/2^+l 


Where  we  define 


=  A+-A_, 


A±=S-lC'lA4±Cj  and 


a±  =  ,±  =  S~XCAXDA±CAS 


The  “hats”  and  “bars”  indicate  that  matrices  use,  respectively,  either  “Roe”  averaged  data  or 
arithmetic  averaged  data.  Note  that  after  some  algebraic  manipulation  the  Modified  Steger 
Warming  method  above  can  be  written  in  a  similar  form  as  the  Roe  method 


or 


1  i+l/2  +/+l/2  /  "*/+l/2  *+l 

F^t"  - (F,„-F,)+  '-a.„n(Fn>  -F,) 

1  7= 


+  2  &+MI2  (Fi+ 1  +  Fi)~  ~^~M/2  (^+1  +  Fi) 

/•«?'-'>  =(SV„„  -a-Mt  )(F„,-F,) 


w 


1+1/2 


p(M-S-W)  _ 1^ 

M/2  2  2 


The  TVD  expression  for  the  second  order  accurate  flux  approximation  using  the  Modified 
Steger-Warming  method  is 


F" 

1  i+l/2 


=  #♦. 


i+l/2 


I-A  — 
1+1/2  Ax 


W-K»+*- 


i+l/2 


I +  A- 


/+ 1/2 


The  matrix  !3l±  can  be  diagonalized  as  follows  S±  -  S  ' Ca'D±CaS  .  Therefore  we  can  write  flux 
vector  as 
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pn  —  CD  Fn  a-  Of  Fn  _i_  J_  c -1  r'~\  n  r* 

-*,+ 1/2  +/+l/2  " /  '  i  '  -  ^  j 


V+l/2  /+!  ^/+l/2^/f,/+l/2A^+,/+l/2'^/l,/+l/2‘J/+l/2 


__  e->  r-1  n  r  s 

2  ‘J/+l/2'^/4,/+l/2  -,/+l/2^/4,/+l/2°/+l/2 


/  1  A' 

I  ~  A*,..-  - 


V 


/+1/2 


Ax 

A/ 


/^,+1/2-l^-C) 


or 


/7W  —  CD  Fn  4-  ^  Z7'7  i _ _  e»  — 1  F) 

^z+1/2  ^  +/+|/2//  "r  ^  ”7+1/2  <+l  “r  2  °'+l/2^/4,/+l/2X>+,/+l/2 


Q(+1/A+,/2(^-^,) 

C^nS^nWi-K) 

Let  the  eigenvalues  of  A  be  ordered  as  A^=u ,  A^=u  +c  and  =  u  -c  ,  then 


_  J_  c-'  r-'  n 

2  0i+l/2'^/f,i+l/2i-y-,;+l/2 


( 

A/  ^ 

I-\ 

V 

^ +,/+!/ 2 1 

AX  , 

/ 

1  ^ 

I- 

V 

1^-,  z+l/2 

1  AX, 

X. 

0 

0  ' 

_  1 

~2 

l±sgn(4) 

0 

0 

D±  = 

0 

0 

<2 

0 

0 

<3. 

0 

0 

l±sgn(^) 

0 

0 

l±sgn(^)_ 

We  define  the  the  vector  Gk  for  k  =  i-\, i, i  + 1  and  i  +  2  to/  =  k  + 1 ,  by 


Gk  = 


S  i 

Si 

Si. 


—  C  C  Fn 
~  A,+in  Mn  k  ’ 


Then 


1  - 


Kt/i  -  S,'\hGa.mh  {^+.,+1/2  G,  +  D.j+U2G,  +  -  D+  m/2 


- D 


-.1+1/2 


( 

A/" 

/- 

V 

^z+1/2  | 

Ax, 

f 

1  ^ 

I- 

V 

|“^  z+l/2 

1  Ax, 

(G-G„) 
&m-Gm)  } 


We  now  apply  the  flux  limiter  to  define  the  TVD  flux  for  the  Modified  Steger-Warming  method. 

A/ 


Fn  -  c->  c 

1  /+ 1/2  “  ° i+l/ 2^ 4 


z+l/2 


d+.iSij+d-*guM  +  1  - Kj+i/2 |— J(^+.|V'+  (s-!., -&./-i)  -d-iVteu+\  -£|,)) 
d+,iSl.i  +  d-,\Sl,i+\  +  2^M  —  I /^2,/+l/2 1  Av  )(^/+,2,//+(^2,;  —  S 2, /-l )  _  2*^+  (^2,1+1  —  &2./)) 


A/ 
'Ax. 

I At} 
'  Ax 


+  d  1-|4w/2|t^-  \(d+.W+(Si,t  -  Si,^)- d-^Si,M-  SXi)) 


where 


wASij-Sij- i)  =  min  m°d(g/., -Si,i-i>S,,M -Si.,) 
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An  example  of  a  program  for  calculating  the  flux  F"+Xll  at  all  interior  surfaces,  from  /'  =  1  +  1/2  to 
/  =  /  — 1/2,  using  temporary  variables  h \, ,  for  /  =  1, 2  and  3 ,  and  vectors  Gk ,  with  elements  gk , , 
for  k  =  1, 2, 3  and  4  and  /  =  1, 2  and  3 ,  etc,  follows. 


1)  Begin  for  /  =  1, •••,/- 1  do  the  following 

2)  Calculate  G,  =  CA  SMnF”_k  ,  where  i'  =  max(2,i) 

3) Calculat  eG^y,/ 

4)  Calculate  G^C^S^F^ 

5)  Calculate  G4  =  CA  v2Si+U2F"+  2,  where  /"  =  min(/-2,/) 


6)  Calculate,  for  /  =  1, 2  and  3 , 


if  sgn(i,/+1/2)  =  1  then  h,  =  gv  +  ^f  1-|^,/+1/2|^- )minmod(g2/ -gu„g3J -g2J) 


AxJ 
A/a 
Ax 


else  h,=gy-U  l-|4,i+1/2|A  minmod(g3/-gw,g4,-g3/)  , 


7)  Calculate  Ff+Xll= 


K 

K 

K 


8)  End  if  i  =  I-l 


The  figure  below  shows  the  results  for  the  Modified  Steger-Warming  method  with  TVD  for  flow 
within  a  shock  tube. 


Figure  6  Density  comparison  for  the  Modified  Steger-Warming  method,  with  TVD,  for  a 
moving  shock  wave,  contact  discontinuity  and  rarefaction,  (CFL=0.9,  40  time  steps) 
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III.  Applications 

(1)  External  Flow:  Simulation  of  the  Ziemer  Experiment 


In  1959  R.W.  Ziemef'  reported  results  from 
an  experimental  investigation  in  magneto¬ 
aerodynamics.  He  placed  a  hemi-spherically 
nosed  cylinder,  of  diameter  0.02m  and  made 
of  Pyrex  glass,  within  an  electromagnetic 
shock  tube  producing  a  hypervelocity  flow 
of  ionized  air.  He  observed  that  with  the 
magnetic  field  turned  on  the  shock  wave 
standoff  distance  increased  by  a  factor  of  7.5 
for  a  magnetic  interaction  parameter 
Q  =  <J\Blrh(ly I pmum  =  69  .The  assumed  free 
stream  conditions  are 


velocity 

5690  m/s 

pressure 

3033  NW 

temperature 

9813  °K 

density 

5.847x1 0'^kgW 

Figure  7  Pressure  contours  within  magnetic 
field  Bo=2T 


A  discharge  of  electric  current  was  used  first  to  ionize  the  air  at  one  end  of  a  3  inch  diameter 
Pyrex  glass  pipe  and  was  then  followed  by  a  further  large  capacitor  discharge  through  the 
ionized  gas.  The  large  joule  heating  of  the  gas,  plus  the  self  pinching  of  the  gas  by  the  large 
magnetic  field  created  by  the  current  flow,  produced  a  strong  shock  wave  that  traveled  down  the 
glass  pipe  toward  the  test  section.  A  high  speed  camera  was  used  to  record  the  passage  of  the 
shock  front  past  the  model  and  during  the  establishment  of  steady  equilibrium  gas  flow.  A  bow 
shock  wave  formed,  dividing  the  shock  layer  about  the  model  from  free  stream  ahead.  A  dipole¬ 
like  magnetic  field  was  created  by  a  current  pulse  through  the  copper  coil  located  within  the 
model,  which  interacted  with  the  surrounding  ionized  flow. 

The  magnetic  interaction  parameter  is  defined  by  Q  =  oxBlrbdyl pjim,  where  cr,  is  the 
electrical  conductivity  of  the  gas  within  the  shock  layer,  B0  is  the  strength  of  the  dipole  at  the 
stagnation  point,  rhJy  is  the  body  radius,  pn  is  the  free  stream  density  and  u m  is  the  free  stream 

velocity.  The  shock  wave  standoff  distance  was  observed  to  increase  by  a  factor  of  7.5  with 
Q  =69  over  that  with  B0  =  0  . 

Bush  Analysis 

Bushh  cleverly  simplified  the  set  of  governing  equations  into  an  ordinary  differential  equation 
and  found  a  similarity  solution.  He  modeled  the  low  subsonic  flow  in  the  stagnation  region  of 
the  shock  layer  as  incompressible,  assumed  a  spherical  shock  wave  shape,  initialized  the  solution 
at  the  shock  and  marched  it  toward  the  body.  He  also  assumed  that  the  electrical  conductivity 
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was  constant  within  the  shock  layer  and  zero  outside  of  it.  His  analytical  results  supported 
Ziemer’s  experimental  observations  and  predicted  a  reduction  in  heat  transfer  through  magnetic 
field  interaction. 

Poggie  and  Gaitonde  Computation 

Poggie  and  Gaitonde7  used  the  theoretical  analysis  of  Bush  for  validating  their  computer 
program  developed  for  solving  the  equations  of  magneto-fluid  dynamics  using  the  low  magnetic 
Reynolds  number  approach.  This  approach  fixes  the  magnetic  field  to  that  imposed.  The 
magnetic  field  interacts  with  the  flow  through  Lorentz  forces  and  joule  heating,  modeled  through 
source  terms  added  to  the  Navier-Stokes  equations.  Their  simulations,  using  perfect  gas 
relations,  agreed  well  with  Bush’s  theory,  thus  providing  validation  for  their  numerical 
procedures. 

Present  Study 

The  purpose  of  the  present  study  is  the  development  of  numerical  procedures  to  simulate 
the  interaction  of  strong  magnetic  and  electric  fields  with  weakly  ionized  aerodynamic  flows. 
The  current  within  the  fluid  can  induce  an  additional  magnetic  field,  which  is  governed  by  the  set 
of  Maxwell  equations.  The  complete  Navier-Stokes  equations,  with  magnetic  stress  and  energy 
terms  included,  plus  Maxwell’s  equations  are  solved  together,  subject  to  the  usual  MFD 
assumptions  of  charge  neutrality  and  that  the  current  density  is  proportional  to  the  curl  of  the 
magnetic  field.  This  set  will  be  called  the  full  equations  of  magneto-fluid  dynamics  (MFD) 
herein,  although  a  larger  set  without  the  MHD  approximations  is  more  accurately  termed  the  full 
set.  Ziemer’s  experiment  is  used  herein  as  a  difficult  challenge  for  algorithm  development  and 
flow  simulation.  Although  the  simulation  to  be  presented  is  expected  to  be  closer  to  reality  than 
those  of  Bush  and  Poggie  and  Gaitonde,  it  is  still  limited  as  will  be  explained  below. 

The  full  set  of  MFD  equations,  Ohm’s  for  conducting  gases,  and  the  equilibrium  chemistry 
model  were  used  in  the  following  simulations.  A  blast  wave  moved  at  Mach  21.5  into  stationary 
air,  at  temperature  273°  K  and  pressure  9.33N/m2,  past  the  model.  The  free  stream  Mach  number 
for  the  flow  behind  the  blast  wave,  relative  to  the  model,  was  of  2.0.  The  free  stream  conditions 
are  given  in  the  table  below.  Although  an  attempt  was  made  to  match  the  experimental 
conditions  of  Ziemer,  no  precise  match  for  the  magnetic  interaction  parameter  Q  could  be  made 
because  of  the  large  variation  in  electrical  conductivity  behind  the  shock.  There  was  no  clear 
choice  in  which  value  to  use.  The  imposed  magnetic  field  at  the  stagnation  point  was  B0= 0,  1  or 
2  Tesla.  The  temperature  behind  the  shock  wave  near  the  nose  was  about  23,000  °K  and 
pressures  were  as  high  as  2xl04  N/m2.  At  these  conditions  the  flow  was  almost  completely 
dissociated.  The  chemistry  model,  shown  above,  should  also  have  included  doubly  ionized 
molecules  and  atoms,  which  limits  to  some  degree  the  realism  of  the  present  simulation.  The 
present  chemistry  model  should  be  sufficient  for  temperatures  less  than  1 5,000  °K,  which  would 
cover  the  main  region  of  current  hypersonic  flow  interest.  Because  the  main  purpose  of  the 
present  simulation  is  to  challenge  the  algorithm  in  development  and  because  the  extension  to 
doubly  ionized  atoms  and  molecules  is  very  challenging  in  itself,  the  chemistry  model  has  not  at 
present  been  further  extended.  An  isothermal  wall  boundary  condition  was  assumed  with  the 
wall  temperature  at  273°K. 
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The  equation  for  the  magnetic  dipole  is 

r3  r3 

given  by  B  =  B0 -^j-cos0er  +  B0 -^-sin  9ee , 
r  2  r 

where  r  is  measured  from  the  origin,  6  is 
measured  from  the  negative  x-axis  direction, 
B0  is  the  dipole  strength  and  rMv  is  the  body 

radius.  The  dipole  should  fall  off  in  strength 
away  from  the  body  along  the  axis  of 

symmetry  as  B  =  B0  [rMy  / r)  .  However, 

Ziemer  measured  a  faster  fall  off  rate  of 

/  \3-61 . 

B  =  B0 1  rMy  /  r  I  in  his  experiment. 


Figure  8  Dipole  magnetic  field  lines 


Stando[[Distance 

The  pressure  contours  for  flow  about  the  model  are  shown  below  in  Figs. 9-1 1,  with  the  results 
for  the  full  set  of  MFD  equations  at  the  top  half  of  each  figure  and  those  for  the  low  magnetic 
Reynolds  number  approach  below.  The  standoff  distance  of  the  bow  shock  wave  from  the  model 
is  compared  for  each  approach  versus  magnetic  dipole  strength  in  Fig.  12,  with  the  dashed  curve 
indicating  the  low  magnetic  Reynolds  number  approach. 


Figure  9  Pressure  contours,  B0= 0 


Figure  10  Pressure  contours,  50=1 


Assuming  that  the  numerical  simulations  are  realistic,  it  is  clear  from  the  above  figures  that 
the  two  approaches  yield  slightly  different  standoff  distances  and  pressure  contour  results,  as 
seen  in  Figs.  10  and  11.  To  explain  this  difference  we  can  observe  the  imposed  and  induced 
magnetic  fields  shown  below  in  Figs. 13-16.  The  maximum  induced  magnetic  field  is  less  than 
5%  of  the  imposed  magnetic  field  dipole  strength  and  it  falls  off  at  a  slower  rate  with  distance 
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from  the  body  The  induced  field  is  in  general  of  opposite  sign  to  that  imposed,  which  in  effect 
cancels  part  the  imposed  field  and  lessens  its  effect  on  the  bow  shock  wave  standoff  distance. 


x 

Figure  11  Pressure  contours,  B0= 2 
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Figure  13  Imposed  Bx  field,  B0= 2 


Figure  15  Imposed  By  field,  B0= 2 
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Figure  14  Induced  Bx  field,  B0= 2 
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Figure  16  Induced  B  field,  B0= 2 
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Fig.  17  shows  the  induced  magnetic  field  lines  for  the  imposed  dipole  strength  B0=  1  case.  Note 
the  similarity  with  the  streamlines  shown  in  Fig.  18.  The  dark  lines  in  these  two  figures  are  Mach 
contours  used  to  indicate  the  position  of  the  bow  shock  wave  and  regions  of  low  speed  flow.  The 
Mach  1.91  is  used  to  indicate  the  shock  wave,  but  it  also  renters  the  figure  downstream  of  the 
shock  wave  as  a  simple  Mach  contour.  The  Mach  0.1  contour  appears  in  Fig.  18.  Streamlines 
entering  within  it  are  not  accurately  placed. 


Compare  the  direction  of  the  magnetic  field  lines  in  Fig.  17  for  the  induced  field  with  those  of 
Fig. 8  for  the  imposed  field.  Note  they  are  opposite  in  direction  ahead  of  the  body.  The  Faraday 
equation  can  be  written  as 


SB  - 
—  =  Vx 
St 


«x5 - ! — Vx5 

production 


diffusion 


(4) 


The  production  term  above  will  continue  to  change  the  induced  magnetic  field  until  m  x  5  =  6  , 
that  is  until  the  flow  and  combined  magnetic  field  are  aligned.  Note  the  near  alignment  between 
the  two  in  Figs.  17  and  18.  The  pinching  of  the  induced  magnetic  field  lines  in  Fig.  17  counteracts 
the  expansion  of  those  in  Fig. 8.  The  diffusion  term  in  Eq.(4)  is  exceedingly  strong  for  most 
aerodynamic  flows  because  of  the  in  general  small  values  of  the  electrical  conductivity. 

The.  Hall  Sli£_  Effect 

The  Hall  and  ion  slip  effects  tend  to  reduce  the  scalar  electrical  conductivity,  as  seen  in 
Eq.(l).  This  can  significantly  reduce  the  interaction  of  the  magnetic  field  with  the  flow.  Figs.  19 
and  20  show  the  pressure  contours  for  B0= 2  for  both  the  full  MFD  and  low  magnetic  Re 
approaches,  with  the  Hall  and  ion  slip  results  at  the  top  half  of  the  figures  and  the  scalar 
electrical  conductivity  results  at  the  bottom  half  of  the  figures. 
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Figure  19  Pressure  contours  with  and  without 
Hall  and  ion  slip  effect,  full  MFD,  B0= 2 


Figure  20  Pressure  contours  with  and  without 
Hall  and  ion  slip  effect,  low  mag.  Re,  B0= 2 


Dr  an 


Figure  21  Normal  surface  stress,  full  Figure  22  Normal  surface  stress,  low  mag. 

MFD  Eqs.  Re  Eqs. 


The  surface  stress,  normalized  by  the  dynamic  pressure  of  the  free  stream,  versus  distance 
along  the  body  from  the  stagnation  point,  is  shown  in  Fig.2 1  for  the  full  set  of  MFD  equations. 
The  hemisphere-cylinder  junction  is  at  s=0.0157.  These  equations  are  written  in  “conservation 
law  form”  (divergence  form)  and  the  force  acting  on  the  body  surface,  both  aerodynamic 
(pressure  and  viscous  stress)  and  magnetic,  are  applied  directly  at  the  body  surface.  On  the  other 
hand,  the  set  of  equations  for  the  low  magnetic  Reynolds  number  approach,  using  source  terms 
to  include  the  magnetic  effects,  applies  magnetic  forces  on  the  body  as  “acting  at  a  distance”. 
Fig. 22  shows  the  normal  stress,  along  the  body  surface  for  the  low  magnetic  Reynolds  number 
approach.  Viscous  stress  was  not  significant  for  these  simulations.  Notice  that  the  magnitudes  of 
normal  surface  stress  for  the  two  approaches  are  widely  different.  One  could  misinterpret  drag 

26 


acting  on  the  body  from  consideration  of  these  two  figures  only.  The  full  MFD  equation 
approach  correctly  predicts  increased  drag  with  increased  imposed  magnetic  dipole  strength. 
Fig.22  for  the  low  magnetic  Reynolds  number  approach  does  not  show  this  same  trend  because 
most  of  the  magnetic  force  drag  occurs  via  “action  at  a  distance”.  Note  the  large  change  in  scale 
of  the  stress  axes  in  these  two  figures  below. 

The  drag  coefficients  for  each  approach  were  calculated  as  follows. 

For  the  full  MFD  equation  approach 


total  axial  surface  stress 


“ MFD 


o  ulxrlfy 


and 


For  the  low  magnetic  Reynolds  number  approach 


11  low  Re 


mag 


total  axial  surface  stress  +  total  axial  volume  force 


where  the  “total  axial  volume  force”  includes  the  Lorentz  source  terms  integrated  over  the  fluid 
volume  surrounding  the  body.  The  Lorentz  force  field  is  shown  near  the  stagnation  point  in 
Fig.23  for  the  low  magnetic  Reynolds  number  approach.  It  appears  to  focus  the  flow  toward  the 
stagnation  point  and  trap  it  ahead  of  the  body.  The  c(,  results  are  shown  for  both  approaches  in 

Fig. 24,  with  the  low  magnetic  Reynolds  number  approach  results  connected  by  the  dashed  line. 
The  figure  shows  that  the  drag  coefficient  is  proportional  to  the  increase  in  the  shock  layer 
volume,  caused  by  the  increase  in  shock  standoff  distance. 


Figure  23  Lorentz  force  field,  low  Remag  Eqs.,  B0=l  Figure  24  cd  vs.  dipole  strength 
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Heal  Transfer 

As  can  be  seen  in  Figs. 25  and  26,  the  two  approaches  predict  different  rates  of  heat  transfer, 
though  both  predict  decreases  in  heat  transferred  to  the  body  with  increased  imposed  magnetic 
dipole  strength.  The  normalized  values  of  heat  transfer  plotted  in  the  figures  are  given  by 

k  — 

2  Ao«« 


where  n  is  normal  to  the  surface  coordinate.  The  difference  in  heat  transfer  for  dipole  strengths 
B0=l  and  B0=2  appear  smaller  than  that  for  either  one  and  Bo=0.  The  stagnation  point  heat 
transfer  is  significantly  reduced  for  the  full  MFD  approach,  but  not  so  for  low  magnetic 
Reynolds  number  approach,  which  is  at  variance  with  the  analysis  of  Bush  and  the  simulations  of 
Poggie  and  Gaitonde.  This  remains  unexplained  at  present. 


Figure  18  Surface  heat  transfer,  full  Figure  19  Surface  heat  transfer,  low 

MFD  Eqs.  Remaf  Eqs. 


Free  Stream  Effects 

It  is  thought  that  the  effect  of  the  Lorentz  force  field  upon  the  free  stream,  ahead  of  the  bow 
shock  wave,  may  important.  The  work  done  on  the  fluid  by  these  forces  plus  joule  heating  could 
add  heat  to  the  free  stream,  consequently  raising  its  temperature  and  increasing  the  speed  of 
sound.  This  effectively  would  lower  the  free  stream  Mach  number,  and  perhaps  explain  the 
increase  in  standoff  distance.  The  Bush  analysis  assumed  that  the  electrical  conductivity  outside 
the  shock  layer  was  zero.  Fig.27  compares  the  flow  fields:  (1)  with  the  electrical  conductivity  set 
to  zero  ahead  of  the  bow  shock  wave,  top  half,  and  (2)  with  it  defined  as  before,  bottom  half. 
Both  used  the  low  magnetic  Reynolds  number  approach  with  B0=2  Tesla.  There  is  little 
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difference  in  the  results,  as  shown  in  the  figure,  indicating  that  the  magnetic  field’s  effect  on  the 
free  stream  is  not  important  for  this  case  and  Bush’s  assumption  was  valid. 


Figure  27  Pressure  contours,  low  Rema?  Eqs,  Figure  28  Pressure  contours,  B0= 2,  full 
B0= 2,  <re  =  0  outside  shock  (bottom  half)  MFD  (top  half),  low  Remu>,  (bottom  half) 

Full  MFD  vs.  low  Re 

-  mag 

From  the  above  results  and  discussion,  it  can  be  thought  that  the  difference  between  the  two 
approaches  is  solely  caused  by  the  differences  in  the  magnetic  fields,  caused  by  the  addition  or 
not  of  the  induced  field  to  the  imposed  field.  This  can  be  tested  by  restarting  a  converged  flow 
simulation  using  the  full  MFD  approach  as  the  initial  condition  for  a  low  magnetic  Reynolds 
number  simulation.  If  the  above  supposition  is  true,  the  continued  simulation  will  remain 
unchanged.  Fig.28  shows  the  pressure  contours  for  both  approaches,  the  converged  MFD  results 
on  top  and  the  continued  solution,  using  low  magnetic  Reynolds  number  approach  results,  below. 
Although  the  shock  standoff  distance  remains  the  same,  the  contour  patterns  are  somewhat 
different.  This  was  not  expected  and  remains  unexplained  at  present. 


(2)  Internal  Flow:  Simulation  of  an  “Energy  Bypass”  Scramjet  Engine 


Figure  29  Generic  Hypersonic  Vehicle 


Figure  30  Starboard  Side  Engine  Box 


Fig.29  shows  the  surface  heating  about  a  generic  hypersonic  vehicle  traveling  through  the 
atmosphere  at  Mach  25.  The  heating  is  extreme  at  the  vehicle  nose,  wing  leading  edges  and 
engine  inlet.  MFD  may  perhaps  be  used  to  control  the  flow  in  these  regions  and  possibly 
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alleviate  severe  heating.  It  also  offers  potential  to  thrust  enhancement  for  proposed  scram  jet 
engines.  The  starboard  side  engine  box  of  the  above  vehicle  is  shown  in  Fig.30. 


To  test  the  new  algorithm  features  discussed  herein  (the  Hall  Effect,  ion  slip  and  calculation  of 
electrical  conductivity  via  electron  and  ion  number  densities),  a  simplified  test  geometry  was 
chosen  for  the  engine  box.  It  consists  of  a  rectangular  geometry  with  three  sections  for  MFD 
power  generation,  combustion  and  MFD  acceleration,  as  shown  in  Fig.3 1 .  Each  section  is  3m 
long,  5m  wide  and  lm  deep. 


Magnetic  and  electric  fields  are  imposed  on  the  flow  normal  to  the  flow  direction  and  to  each 
other.  The  imposed  Magnetic  field  was  Bv  =  1  Tesla  and  a  local  loading  factor  of  0.75  was  used 

to  decelerate  the  flow  in  the  generator  section,  from  an  entering  flow  of  Mach  10  to  about  3.5  at 
the  combustor  entrance.  Heat  was  added  within  the  combustor  at  a  rate  proportional  to  a  fraction 
of  the  kinetic  energy  of  the  uniform  shock-heated  free  stream,  assumed  ahead  of  the  inlet, 
(w0  =  6304/w/s,  p0  =  60  N  /  m2  and  T0  =  1000°  K).  The  magnetic  field  imposed  upon  the  flow 


within  the  accelerator  was  Bv  =  0.25  Tesla  and  the  load  factor  based  on  free  stream  velocity  was 
2. 
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Figure  32  Induced  Magnetic  Field  Components 
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The  induced  magnetic  field  component  contours  are  shown  in  Fig.32  on  an  x-y  plane  at  the 
center  of  the  channel. 
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Figure  33  Velocity  Vectors  and  u-Velocity  Contours 


Velocity  vectors  and  contours  of  the  axial  velocity  component  are  shown  in  Fig.33.  The 
deceleration  of  the  flow  within  the  MFD  generator  section  is  shown,  followed  by  further 
reduction  within  the  combustor  and  then  acceleration  within  the  MFD  accelerator  section. 


Contours  for  electron  concentrations  and  electrical  conductivity  are  shown  in  Figs.  34  and  35. 
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Figure  34  Electron  Concentration  Contours 


Figure  35  Electrical  Conductivity  Contours 


Figure  36  Specie  Concentrations  Figure  37  Electron  Concentration 
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I 


► 


Specie  concentrations  for  N2,  02,  NO,  N,  O  and  e  half  way  through  the  accelerator  section, 

at  x=7.5m  ,  are  shown  in  Figs.  36  and  37.  Fig. 38  shows  that  both  axial  thrust  and  maximum 
velocity  through  the  channel,  normalized  by  values  at  the  entrance,  both  decrease  with  power 
generation  and  increase  by  almost  a  factor  of  50%  within  the  accelerator. 
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Figure  38  Axial  Thrust  and  Maximum  Velocity  vs.  x 


Conclusion 

Algorithms  have  been  developed  to  simulate  weakly  ionized  flows  relevant  to  flow  about 
hypersonic  flight  vehicles  and  through  energy  bypass  scram  jet  engines.  These  include  numerical 
procedures  for  casting  the  governing  MFD  equations  in  a  form  enhancing  numerical  precision, 
boundary  condition  application,  scalar  and  tensor  formulations  for  electrical  conductivities  for 
Hall  current  and  ion  slip  effects,  calculation  of  equilibrium  properties  for  an  open  ended  set 
component  gas  species,  control  of  magnetic  field  divergence  and  improving  the  fully  implicit 
algorithm  for  solving  the  governing  equations.  Fairly  realistic  applications  were  made  test  the 
robustness  of  the  developed  procedures. 

Although  the  main  thrust  of  the  present  study  is  algorithm  development  for  weakly  ionized 
hypersonic  flow  simulations  within  strong  magnetic  and  electric  fields,  fairly  realistic 
simulations  were  made  of  the  Ziemer  experiment,  conducted  in  the  late  1950s.  This  is  a  good  test 
problem  for  algorithm  development  because  of  the  strong  interaction  of  the  magnetic  field  with 
the  flow  and  because  the  gas  temperatures  are  high  enough  to  simulate  hypersonic  conditions. 
The  present  results  confirmed  earlier  studies  predicting  increased  drag  and  decreased  heat 
transfer  with  increased  magnetic  field  interaction.  They  also  explain  why  the  low  magnetic 
Reynolds  number  approach  to  flow  simulation  predicts  larger  bow  shock  wave  standoff  distances 
than  the  full  MFD  equation  approach,  because  the  induced  magnetic  field  tends  to  partly  cancel 
the  imposed  magnetic  field.  In  general  the  two  approaches  agree  well  with  one  another. 
Unexplained  differences  remain  between  the  two  approaches  that  invite  further  study.  The  Hall 
and  ion  slip  effects  significantly  reduce  the  interaction  of  the  magnetic  and  flow  fields. 
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The  numerical  procedures  were  applied  to  simulate  the  flow  through  an  energy  bypass  scram  jet 
engine.  The  gas  was  assumed  to  be  the  chemical  equilibrium  containing  1 3  species  and  1 1 
chemical  reactions.  The  electron  and  ion  concentrations  were  the  used  to  determine  electrical 
conductivities  of  seeded  air.  The  algorithm  has  been  applied  to  a  simplified  rectangular  channel 
geometry  containing  sections  for  MFD  power  generation,  combustion  and  MFD  acceleration. 
Predictions  of  engine  thrust  were  calculated. 
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